/*
  This function calculates one-dimensional interpolation with "norder"-th order accuracy.
  i.e. if norder = 1 : linear interpolation etc.

  It performs odd-order interpolation only. i.e. if norder = 2, 4, 6 etc., error will be returned.

  It accepts "nxpt" number of grid-points "x", and "nxiplpt" number of interpolation-points "xiplpt" (where the
  one-dimensional interpolation is to be performed. It returns the co-efficient matrix "mtrx" of size [nxiplpt]x[nxpt].
  It doesn't filter out the zero-elements. */

# include <stdlib.h>
# include <math.h>

static void lagrange_interpol(double *ipval, double *xstn, unsigned short nstn, short ipstn, double dxstn);

/* This evaluates the co-efficient matrix for norder-th order interpolation in one-dimension. */
extern short cmtrx_interpol(double *x, unsigned short nxpt, double *xipl, unsigned short nxiplpt,\
 unsigned short norder, double **mtrx)
  {
   unsigned short j, *tmtrx;
   short k, idxmin;
   double dx, dxmin, epsilon;

/* Return error if order of the interpolation is even. */
   if(norder%2 == 0) return -10;

/* This the transformation matrix and it is used to mark the beginning of the stencil
   on which the interpolation will be carried out. It becomes very useful when we
   approach towards ends of the domain. */
   tmtrx = (unsigned short *)malloc(sizeof(unsigned short)*nxiplpt);

/* This is used to check if two float quantities are same or not. */
   for(epsilon = fabs(x[nxpt-1]-x[0]), j = 1; j < nxpt; j++) if(epsilon > (dx = fabs(x[j]-x[j-1])) ) epsilon = dx;
   epsilon /= 100;

   for(k = 0; k < nxiplpt; k++)
      {
    /* Finding out the closest grid-point from the k-th interpolation point. */
       dxmin = x[nxpt-1]-x[0];
       idxmin = -1;
       for(j = 0; j < nxpt; j++)
          {
           mtrx[k][j] = (double)0; /* Initializing the matrix's k-th row with zero. */

           if(fabs(dx = xipl[k]-x[j]) < fabs(dxmin))
              {/* Locating the closest grid-point from the k-th interpolation-point. */
               dxmin = dx; /* The closest distance. */
               idxmin = j; /* The grid-number corresponding to the closest grid. */
              }
           }

   /* Check if the k-th interpolation-point coincides with any of the grid-points. */
      if(fabs(dxmin) > epsilon)
          {/* The k-th interpolation-point doesn't coincide with a grid-point. */
           if(nxpt <= norder) return -11; //Return error if number of points are insufficient.

        /* Get the grid-number right before the k-th interpolation point. */
           if(dxmin < 0) idxmin--;

        /* Mark the beginning of the stencil that will be used in interpolation. */
        /* The k-th interpolation point is close to beginning of the grids. */
           if(idxmin <= (norder-1)/2) tmtrx[k] = 0;
        /* The k-th interpolation point is close to end of the grids. */
           else if(idxmin >= nxpt-(norder+1)/2-1) tmtrx[k] = nxpt-norder-1;
        /* The k-th interpolation point lies in middle of the grids. */
           else tmtrx[k] = idxmin - (norder+1)/2 + 1;

        /* This interpolates using the Lagrange-polynomial at the k-th interpolation-point. */
           lagrange_interpol(&mtrx[k][tmtrx[k]], &x[tmtrx[k]], norder+1, -1, xipl[k]-x[tmtrx[k]]);
          }
   /* The k-th interpolation-point coincides with a grid-point. */
      else mtrx[k][idxmin] = (double)1;
     }

   free(tmtrx);

   return 10;
  }



/* This function calculates a row-vector of the co-efficient matrix responsible for evaluation of the
  one-dimensional interpolation at a given point. It usage the Lagrange-polynomial approach and th
  row-vector consists of the Lagrange-weights at each stancil-point.

  xstn: Beginning of the stencil-points.
  nstn: Number of stencil-points. i.e. stencil-size.
  ipstn: If the point at which the interpolation is done coincides with grid-points,
         it marks that grid-point on the stencil. "dxstn" is ignored in this case.
  dxstn: If the interpolation-point doesn't coincide with grid-points, it marks the distance
         of the point on the stencil. ipstn = -1.

  "ipstn" and "dxstn", both are relative to beginning of the stancil. */

static void lagrange_interpol(double *iplval, double *xstn, unsigned short nstn, short ipstn, double dxstn)
  {
   unsigned short i, j;

   if(ipstn < 0)
      {
       for(j = 0; j < nstn; j++)
         for(iplval[j] = (double)1, i = 0; i < nstn; i++)
           if(i != j) iplval[j] *= (xstn[0]+dxstn-xstn[i])/(xstn[j]-xstn[i]);
      }
   else for(j = 0; j < nstn; j++) iplval[j] = (double)(j == ipstn ? 1:0);

   return;
  }
